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Abstract. The combinatorial explosion of empirical parameters in tens of thousands 
presents a tremendous challenge for extended Debye-Hiickel models to calculate activity 
coefficients of aqueous mixtures of most important salts in chemistry. The explosion of pa¬ 
rameters originates from the phenomenological extension of the Debye-Hiickel theory that 
does not take steric and correlation effects of ions and water into account. In contrast, the 
Poisson-Fermi theory developed in recent years treats ions and water molecules as nonuni¬ 
form hard spheres of any size with interstitial voids and includes ion-water and ion-ion 
correlations. We present a Poisson-Fermi model and numerical methods for calculating the 
individual or mean activity coefficient of electrolyte solutions with any arbitrary number of 
ionic species in a large range of salt concentrations and temperatures. For each activity- 
concentration curve, we show that the Poisson-Fermi model requires only three unchanging 
parameters at most to well fit the corresponding experimental data. The three parameters 
are associated with the Born radius of the solvation energy of an ion in electrolyte solution 
that changes with salt concentrations in a highly nonlinear manner. 


1. INTRODUCTION 


Thermodynamic modeling of ac 
chemical and biological sciences [l 


ueous electrolyte solutions plays an important role in 


13]. Despite intense efforts in the past century, robust 


thermodynamic modeling of electrolyte solutions still presents a difficult challenge and re- 
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mains a remote ambition in the extended Debye-Hiickel (DH) models due to the enormous 
number of parameters that need to be adjusted, carefully and often subjectively HQ- 
For example, the Pitzer model requires 8 parameters for a ternary system and up to 8 tem¬ 
perature coefficient s (p arameters) for every Pitzer parameter in a temperature interval from 
0 to about 200 °C jll JQ. It is indeed a frustrating despair (frustration on p. 11 in |9] and 




jye-Hiickcl theory Q] 
5] although potentials 


despair on p. 301 in lj) that approximately 22,000 parameters for combinatorial solutions 
of the most important 28 cations and 16 anions in salt chemistry have to be extracted from 
the available experimental data for one temperature ll|. The Pitzer model is still the most 
widely used DH model with unmatched precision for modeling aqueous electrolyte solutions 
over wide ranges of composition, temperature, and pressure (l3|. 

The Pitzer model and its variants 13] are all derived from the De 
that in turn is based on a linear Poisson-Boltzmann (PB) equation 
calculated from PB near ions (for example) are often far beyond the linear range of the po¬ 
tential near ions or interfaces. The PB equation treats ions as point charges without steric 
volumes and water molecules as a homogeneous dielectric medium without steric volumes 
either and with a constant dielectric constant that neglects ion-water and ion-ion correla¬ 
tions. These simplifications give rise to the elegant, simple, and useful DH theory. However, 
it is precisely because of the linearization and simplifications on steric and correlation effects 
that extended DH models have needed an explosion in the number of parameters in order to 
overcome the deficiencies (simplifications) of the classical Poisson-Boltzmann theory. The 
nonlinear PB equation was developed by Gouy and Chapman 151, Il6 |. 

In the past few years, we have intensively investigated these two effects in a range of areas 


from electric double layers 


HQ, 


ion activities 


19], to biological ion channels 


18 


20 


24] 


and consequently developed an advanced theory — the Poisson-Fermi (PF) theory — that 
treats ions and water molecules as nonuniform hard spheres of any size with interstitial voids 
and includes many of the correlation effects of ions and water. We refer to our previous 
papers and references therein for a historical account of the literature of this theory. In 19], 
we proposed a PF model for calculating activity coefficients of individual ions in aqueous 
single NaCl and CaCl 2 electrolyte solutions at the temperature 298.15 K. The model is 
further tested in this paper for eight 1:1 electrolytes (LiCl, LiBr, NaF, NaCl, NaBr, KF, 
KC1, and KBr), six 2:1 electrolytes (MgCl 2 , MgBr 2 , CaCl 2 , CaBr 2 , BaCl 2 , and BaBr 2 ), one 
mixed electrolyte (NaCl + MgCl 2 ), one 1:1 electrolyte (NaCl ) at various temperatures from 
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298.15 to 573.15 K, and one 2:1 electrolyte (MgCl 2 ) at various temperatures from 298.15 to 

523.15 K, for which the experimental data were compiled by Valisko and Boda in 
Rowland et al. in 13] from various experimental sources in 


26 


25 ] and 


-35]. 


The PF model is developed to calculate individual ion activities for which experimental 


measurements and determination 


10 


36 


3.7], interpretation of measurement data 26 


37- 


39], and comparison of different experimental methods [37], |40] have been extensively in¬ 


vestigated by Wilczek-Vera, Rodil, and Vera in the past two decades. PF results on mean 
activity coefficients can be compared with experimental measurements using the Debye- 
Hiickel equation of individual ion activities [fj]. 

In contrast to the Pitzer model, we show that all experimental data sets of individual 
or mean activity coefficients as a function of variable concentration in single electrolytes or 
mixtures at various temperatures can be well fitted by the PF model with only 3 parameters 
at most for each activity-concentration data curve. The model is characterized by three 
different domains, namely, the Born ion, hydration shell, and remaining solvent domains in 
which the Born ion domain is most crucial because all activities around an ion are mainly 
governed by the singular charge of the ion located at the center of the domain. The Born 
ion domain is defined by the Born radius of the solvated ion, which is unknown and changes 
with salt concentrations in a highly nonlinear manner. 

The three parameters characterize three orders of approximation of the Born radius in 
terms of ionic concentrations. Parameter 1 describes a correction of the experimental Born 
radius of a single ion in pure water without any other ions. Parameter 2 describes an 
adjustment of the unknown Born radius in electrolyte solution that accounts for the Debye 
screening effect, which is proportional to the square root of the ionic strength of the solution. 
Parameter 3 is an adjustment in the next order approximation beyond the DH treatment 
of ionic atmosphere. The physical origin of these parameters is clear unlike that of most 
parameters in the Pitzer method [ll, 41]. It may even be possible in later work to calculate 
some of these parameters from more detailed versions of our model. 

Our approach to partition the free energy domain of a solvated ion into the above three 
sub-domains yields a better approximation to calculate the free energy since these sub- 
domains are determined by the experimental data of solvation and thus separate short- and 
long-range interactions of the ion in a more accurate way. This approach nevertheless incurs 
more complicated numerical methods for solving the nonlinear partial differential equations 
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FIG. 1: The model domain 12 is partitioned into the ion domain 0* (with radius Rf orn ), the 
hydration shell domain 12 ,,/* (with radius Rf h ), and the remaining solvent domain 12 s . 


of the PF model in different domains with suitable interface conditions 0 . We therefore 
present numerical methods in detail for future verification and development of the present 
work. 


2. THEORY 


For an aq 


proposed in 


ueous electrolyte solution with K species of ions, the Poisson-Fermi theory 


.181, 21 ] treats all ions and water of any diameter as nonuniform hard spheres 


with interstitial voids between these spheres. The activity coefficient 7 * of an ion of species 
i in the solution describes the deviation of the chemical potential of the ion from ideality 


( 7 i = 1). The excess chemical potential /r| x = /c^Flny* can be calculated by 19, [42] 


A*r = AGi-AG“, AG, = AG“ = I®/(0), 


( 1 ) 


where ks is the Boltzmann constant, T is an absolute temperature, g* is the ionic charge of 
the hydrated ion (also denoted by i), 0(r) is a potential function of spatial variable r in the 
domain 12 = 12* Ul2 s ^ Ul2 s shown in Fig. 1, 12* is the spherical domain occupied by the ion i, 
12 s h is the hydration shell domain of the ion, 12 s is the remaining solvent domain, 0 denotes 
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the center (set to the origin) of the ion, 0(0) is the value of 0(r) at r = 0, and 0°(r) is a 
potential function when the solvent domain does not contain any ions at all with pure 
water only. The potential function 0(r) can be found by solving the Poisson-Fermi equation 

fl 

(4 2 V 2 - 1) V ■ e(r)V0(r) = p(r), (2) 


e r = 


^s ^ ^sh U ^s 


i lc 


6i €j on 6 o in 

p s (r) = Y,k=i QkC k (r) in Q s 
P ( r ) = 0 in Ti sh 

Pi{ r) = qiS( r - 0) in 

C k { r) = Cfc exp (~M{ r) + —S tIC (r) 

V 


2 dj in Q s h U 

0 in 


in Q, 


S (r) = In 


£(r) 

r B 


in Q, 


(3) 

(4) 

(5) 

( 6 ) 


where e 0 is the vacuum permittivity, e w is the dielectric constant of bulk water, e ion is a 
dielectric constant in Qj, a 3 is the radius of a counterion of the ion i, and <5(r — 0) is the 
delta function at the origin. 

The concentration function C k ( r) is described by a Fermi distribution (5), where Cj? is a 
constant bulk concentration for all k = 1, • • • , K + 1, qx+i = 0, f3 k = qk/ksT, v k = 47 ra|/ 3 , 
v 0 = (eE 1 Vk'j /(K +1) an average volume of all kinds of hard spheres, S' trc (r) is called the 
steric potential, T B = 1 — Y^k=i v k^k a constant void fraction, T(r) = 1 — v kC k ( r) 

is a void fraction function, and K +1 denotes water. The radii of fh and the outer boundary 
of Q s h are denoted by Rf orn and Rf\ respectively, whose values will be determined by 
experimental data. It is natural to choose the Born radius Rf orn (not the ionic_radius a^) 


as the radius of fh 


42], We consider both first and second shells of the ion 43, 44], 


The potential 0°(r) (in Eq. (1)) of the ideal system is obtained by setting p s ( r) = 0 in 
(4), i.e., all particles in do not electrostatically interact with each other since q k = 0 for 
all k. The domain is chosen to be sufficiently large so that 0(r) = 0 on the boundary 
of the domain dfl. The ideal potential 0°(r) is then a constant, i.e., A G® is a constant 
reference chemical potential independent of Cf. 

The distribution (5) is of Fermi type since all concentration functions have an upper 
bound, i.e., C k ( r) < l/v k for all particle species with any arbitrary (or even infinite) potential 
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0(r) at any location r in the domain 2lj. The Poisson-Fermi equation (2) and the Fermi 
distribution (5) reduce to the Poisson-Boltzmann equation and the Boltzmann distribution 
when l c = S tTC = 0. i.e., when the correlation and steric effects are not considered. The 
Boltzmann distribution C4(r) = C® exp (— /3k4>(r)) would however diverge if 0(r) tends to 
infinity. This is a major deficiency of PB theory for modeling a system with strong local 
electric fields or interactions 45]. If the correlation length l c ^ 0, the dielectric operator 


?=e.(l-qv 2 ) in Eq. (2) app roximates the permittivity of the bulk solvent and the linear 


response of correlated ions 


17 


20 


46 


47], and yields a dielectric function e(r) as an output 


of solving Eq. (2) [21]. The exact value of F(r) at any r G cannot be obtained from 

Eq. (2) but can be approximated by the simple formula <T(r) f« e;+ CH 2 o( r )(es“ ej)/C§ 2 o 
since the water density function CH 2 o( r ) = Ck+i{ r) is an output of Eq. (5). This formula 


is only for visualizing (approximately) the profile of P or e. 
The input is the correlation length l c in Eq. (3) [l7, 20, 46 


t is not an input of calculation. 


47] . The actual outputs are the 


numerical solutions of the partial differential equations and boundary conditions. 

The factor Vk/v o multiplying the steric potential function S' trc (r) in Eq. (5) is a modi 


21 


hca- 

241 


tion of the unity used in our previous work 19,[21]- The steric energy —^S tlc (r)kBT 
of a type k particle depends not only on the voidness (r(r)) (or equivalently crowding) at r 
but also on the volume Vk of the particle itself. If all Vk are equal (and thus Vk = v 0 ), then 
all particle species at any location r G U have the same steric energy, i.e., uniform 
particles are indistinguishable in steric energy. The steric potential is a mean-held approxi¬ 
mation of Lennard-Jones (L-J) potentials that describe local variations of L-J distances (and 


thus empty voids) between any pair of particles. L-J potentia 


extremely expensive and unstable to compute numerically 


21 


s are highly oscillatory and 
. Calculations that involve 


L-J potentials, or even truncated versions of L-J potentials must be extensively checked to 
be sure that results do not depend on irrelevant parameters. 

3. METHODS 


To avoid large errors in approximation caused by the delta function <5(r — 0) in (4), the 


potential function can be decomposed as HQ 

f d>(r) 

<t>{ r) = 


49] 

0(r) + 0*(r) + 0 L (r) in Oj 


0(r) in fl sh U 


(7) 
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where 0 *(r) = gi/( 47 re; |r — 0 |) and <f>( r) is found by solving 


(l 2 c V 2 - l) V • e s V0(r) = p( r) in U 
—V ■ 6jV0(r) = 0 in hi* 


( 8 ) 

(9) 


without the singular source term pi(r) = g*5(r — 0 ) and with the interface conditions 

for all r G <9h2j, (10) 



0 (r) 

= 0 


e(r)V0(r) • n 


where n is an outward normal unit vector at r G dfli and the jump function [w(r)] = 
lim rs;i _j. r u(r s h) — linir^r u(rj) with r s h G Q s h and r, ; G fli |l3]. The potential function 0 L (r) 
is the solution of the Laplace equation 


V 2 0 L (r) = 0 in 0 , 


( 11 ) 


with the boundary condition 

0 L (r) = 0 *(r) on <9h2j. ( 12 ) 

The evaluation of the Green’s function 0*(r) on dfli always yields hnite numbers and thus 
avoids the singularity in the solution process. The desired solvation energ y A G, in Eq. (1) 
(and thus the individual ionic activity coefficient 7 *) is then evaluated by [171, |49] 


AG, = k B T hi 7 , = -q, 0(0) + ^> L (0) 


(13) 


Since the interface dfli is a sphere centered at the origin, the Laplace potential 0 L (r) = 
qi/{A'Ke i Rf orn ) is a constant in Qj, i.e., Eq. (11) has been exactly solved. 

The Poisson-Fermi equation ( 8 ) is a nonlinear fourth-order partial differential equation 
(PDE) in h2 s . Newton’s iterative method is usually used for solving nonlinear problems. We 


M 


seek a sequence of approximate solutions < 0 m (r) \ by iteratively solving the linearized 

l J m=1 

PF equation 

(l 2 V 2 - l) V • eV0 m - p' s (4>m- 1) 4>m = PsQ>m- 1) - p' s (4>m- 1) 0m- 1 m (14) 

until a tolerable potential function 0 m is reached, where 0 o( r ) is a given initial guess poten¬ 


tial function, p s ( 0 m _ 1 ) = Ef=i QkC™ x (r), C™ : (r) = Cf exp (-&0m-i(r) + ^S'^ c _ 1 (r)), 


vo 


Qtrc _ i„ ( To(r) \ p _ s^K+1 ^ pym-1 n , ( i ^ ( n_ n™-l („\ 
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and p' s ((f)) = J=p s (0). Note that the differentiation in p' s (0) is performed only with respect 
to (f> whereas S trc is treated as another independent variable although S' trc depends on 0 as 
well. Therefore, p'(0) is not exact implying that this is an inexact Newton’s method 50]. 

The fourth-order problem can be resolved by transforming Eq. (14) into two second-order 
PDEs [l7] 

e s (Z 2 V 2 - l) T(r) = p(0 m _ 1 ) in 1 1 sh U 12 a (15) 

CsV 0 m (r) p (0m— i) <f> m (r) e s 4/(r) p (0m— i) dm— i in 12s/i U 12 s (16) 


by introducing a density like variable T = V 2 0 for which the boundary condition is 

T(r) = 0 on <912 s . 


id 


(17) 


Eqs. (9), (15), and (16) are coupled together in the entire domain 12 with the jump conditions 
in (10). Note that linear PDEs (14), (15), and (16) converge to the nonlinear PDE (8) if 
<f>M converges to the exact solution 0 of Eq. (8) as M —» oo, i.e., the approximate potential 
0 m(f) is sufficiently close to the exact potential 0(r) for all r G Q s h U 12 s if the iteration 
number M is sufficiently large (M « 5 to 37 for this work with error tolerance 10^ 3 ). 

The standard 7-point finite difference (FD) method is used to discretize all PDEs (9), 
(15), and (16), where the jump conditions in (10) are handled by the simplified matched 
interface and boundary (SMIB) method proposed in { 17 I ]. For simplicity, the SMIB method 
is illustrated by the following ID linear Poisson equation (in x-axis) 


with the jump condition 


^ = hl ^ 


= — e,—0*(x) at x = £ = <912 i D <912 s , 
ax 


(18) 


(19) 


where 12 = 12jUl2 s , 12,- = (0, £), 12 s = (£, L), f(x) = 0 in 12,, f(x) ^ 0 in 12 s , and 0' = ^0(x). 
The corresponding cases to Eqs. (9), (15), and (16) in y- and z-axis follow in a similar way. 
Let two FD grids points xi and xi+\ across the interface point £ be such that xi < £ < xi + \ 
and £ = (x; + xi + 1)/2 with Ax = xi +1 — xi — 1 A, a uniform mesh, for example, as used in 
this work. The FD equations of the SMIB method at Xi and xi + 1 are 


e; 


7-1 + (2 — Cl) 02 — C20/+1 
Ax 2 

—di4>i + (2 — d2)(f>i+i — 0z+2 

Ax 2 


— fi + 


— fi +i + 


Co 

Ax 2 
do 


Ax 2 ’ 


( 20 ) 

( 21 ) 
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where 


e* - e s 

Ci — —-—, c 2 — 
c* + e s 

2 e 


2e., 


—e,Ax 


d\ = 


Cj + £$ 


, di2 — 


~ f - ^s 


Ci + C s 


Co = 


Ci + c s 
—e.Ax 


do — 


Ci + c s 


0/ is an approximation of 0(xj), and fi = f(xi). Note that the jump value e<j)' at £ is 
calculated exactly since the derivative of 0 * is given analytically. 

Since the steric potential takes particle volumes and voids into account, the shell volume 
V s h of the shell domain Q sh can be determined by Eqs. (5) and ( 6 ) as 

or \ . (v 8 h -v w or 


s% = - In 


v, h c B 


K+l 


= In 


V, h T B 


( 22 ) 


where the occupancy (coordination) number Of is given by experimental data 43,144]- The 
shell radius Rf h of Q s h is thus determined. Note that the shell volume depends not only on 
Of but also on the bulk void fraction 1 B . namely, on all salt and water concentrations 

(Of)- 


As discussed in 25], the solvation free energy of an ion i should vary with salt con¬ 
centrations and can be expressed by a dielectric constant e(Cf ) that depends on the bulk 


concentration Cf of the ion. Therefore, the Born energy 


AG? orn = 


Qt 


(23) 


e w j 87T£ 0 Ri 

with the Born radius Rf in pure water should be modified with the concentration-dependent 
dielectric constant e(Cf). Equivalently, the Born radius in electrolyte solutions can be 
modified from R® by a simple formula 


R 


Born ( 


/_\ 1/2 /_d \ 3/2 

(Cf) = e(Cf )B“, e(Cf)=a\+ai(C i ) +<4 00 


(24) 


where C i = Cf /M is a dimensionless bulk concentration of type i ions, M is the molar 
concentration unit, and aj, a l 2 , and a\ are adjustable parameters for modifying the experi¬ 
mental Born radius R® to fit experimental activity coefficients 7 , that change with the bulk 
concentration conditions Cf of the ion. The Born radii Rf in Table 1 are cited from 25], 
which are computed from the experimental hydration Helmholtz free energies of these ions 
given in {g]. Numerical values in Tables 1 and 2 are all experimental data for which their 
values are kept fixed throughout calculations once chosen. 























Journal of Chemical Physics 
in the press 


10 


The three parameters in Eq. (24) have physical or mathematical meanings unlike many 


parameters in the Pitzer model 41]. Any model or numerical method incurs errors to 
approximate a real system, i.e., it is impossible to obtain real Born radius Rf orn (Cf) exactly. 
The first parameter a\ is an adjustment of the experimental Born radius R { - when Cf = 0 
for all i. The second parameter a l 2 is an adjustment of Rf orn (Cf) that accounts for the real 
thickness of the ionic atmosphere (Debye length), which is proportional to the square root of 
the ionic strength y/1 in the Debye-Hiickel theory jjj] . The third parameter oj s is simply an 
adjustment in the next order approximation beyond the DH treatment of ionic atmosphere. 

We summarize the mathematical solution process for determining the activity of ionic 
solutions in the following algorithm. 


1. Solve Eqs. (9), (10), and (16) for 0 with p' = T = 0 (in pure water), Rf orn = R°, 
and 0 L = q,/(4n6 l R ( -) to obtain AG 9 by Eq. (13) and then set </> 0 = <j>. 

2. Solve Eqs. (15) and (17) for T with Rf orn in (24). 

3. Solve Eqs. (9), (10), and (16) for (f> m with (j) L = qi / (Ane t R,^ orn ) and then set (f) m -i = 4>m- 
Go to 2 until convergence. 

4. Obtain the activity coefficient 7 \ by Eq. (13). 


Table 1. Values of Model Notations 


Symbol 

Meaning 

Value 

Unit 

k b 

Boltzmann constant 

1.38 x 10 " 23 

J/K 

T 

temperature 

Table 2 

K 

e 

proton charge 

1.602 x 10~ 19 

C 

eo 

permittivity of vacuum 

8.85 x 10 ~ 14 

F / cm 

^ion? 

dielectric constants 

1, Table 2 


/ Q - 2 

correlation length 

j = CD etc. 

A 

or 

in Eq. (22) 

18 [43, 44J 


a Li+i a Na+> a K+ 

radii 

0.6, 0.95, 1.33 

A 

a Mg 2 +> a Ca 2 +i a Ba 2 + 

radii 

0.65, 0.99, 1.35 

A 

®F _ 5®C1~ ) a Br“5 a H 2 0 

radii 

1.36, 1.81, 1.95, 1.4 

A 

DO dO dO 

U ‘U+ ’ - rt Na+’ - fX K+ 

Born radii in Eq. (24) 

1.3, 1.618, 1.95 

A 

dO dO 

n Mg 2+1 - ft 'Ca 2 +’ ■ ft Ba 2 + 

Born radii 

1.424, 1.708, 2.03 

A 

C>0 dO dO 

^Cl - ’ 

Born radii 

1.6, 2.266, 2.47 

A 
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C 




c 



FIG. 2: Indivivual activity coefficients of 
mental data [26] oni = Pos + (cation) and 


from 0 to 1.6 M. 


1:1 electrolytes. Comparison of PF results with experi- 
Neg - (anion) activity coefficients 7 , in various [PosNeg] 


Table 2. Values of e w at various T [51] , 
T/K 298.15 373.15 423.15 473.15 523.15 573.15 
e w 78.41 55.51 44.04 38.23 32.23 25.07 


4. RESULTS 


The PF results of ionic activity coefficients for eight 1:1 electrolytes, six 2:1 electrolytes, 
one mixed electrolyte, one 1:1 electrolyte at various temperatures, and one 2:1 electrolyte 


at various temperatures agree with the experimental data 


26 


35] as shown in Figs. 2, 3, 4, 


5, and 6, respectively. The empirical parameters used to fit the experimental data are a\, 
a 2 , and a\ in Eq. (24), whose values are given in Table 3 from which we observe that the 
PF model requires only one to three parameters to fit those data. 

The mean activity coefficient 'fposNeg of a salt Pos p Neg g is calculated via the formula 
In 7 PosNeg = ^ In 7 Pos + ^ In 7 Neg H , where 7 Pos and 7 Neg are individual activity coeffi¬ 
cients obtained by Eq. (13) for each i = Pos and Neg. For the mean activity coefficients 
of either ternary (Fig. 4) or binary (Figs. 5 and 6) systems, we only need to adjust 3 
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0 0.5 1 1.50 0.5 1 1.5 

([PosNeg 2 ]/M ) 1/2 ([PosNeg 2 ]/M ) 1/2 


FIG. 3: Indivivual 
mental data |26| on 


from 0 to 1.5 M. 


activity coefficients of 2:1 electrolytes. Comparison of PF 
i = Pos 2+ (cation) and Neg - (anion) activity coefficients 7 i 


results with experi- 
in various [PosNeg 2 ] 


parameters of one cation (not all ions) as shown in Table 3. 

The activity coefficients by the PF model are quite successful over a large range of tem¬ 
peratures and concentrations as shown in Figs. 4-6. We used the code of the density model 


developed by Mao and Duan 52] to convert the concentration unit from molality (mol. 
kg -1 ) to molarity (M = mol. dm” 3 ) by the standard formula as given in 52| . where the 
density model has been compared with thousands of measurements at high accuracy. The 
pressure values needed in the code at the corresponding temperatures were set to P = (A) 
1.01 (B) 1.01 (C) 15.48 (D) 39.59 (E) 80.50 bar for Fig. 4 and (A) 1.01 (B) 1.01 (C) 4.73 


(D) 39.50 bar for Fig. 5. In Fig. 6 , the ionic strength / = JT Cfzf and the ionic strength 
fraction 7/ Mg ci 2 = 3m M gCi 2 /(3m M gci 2 + "iNaCi) with m Mg ci 2 and m NaC i being the molalities of 
MgCl 2 and NaCl in the mixture, respectively, where Zi is the valence of type i ions. 

We observe from Table 3 that the approximate Rf orn (Cf) (with salts) deviates from R 3 
(without salts) only in the second to fourth decimal place, i.e., numerical values of 7 * are 
very sensitive to the decimal order of a\, a 2 , and a\ because the Born radius Rf orn (Cf) is 
very close to the origin 0 at which the singular charge in Pi(r) = qi5{ r — 0) is infinite. The 
approximation of the shell radius Rf h (or the coordination number OJ in Eq. (22)), on the 
other hand, is much less significant than that of Rf orn because the electric potential 0 PF (r) 
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FIG. 4: Mean activity coefficients of mixed electrolytes. Comparison of PF results (curve) with 


experimental data (symbols) compiled in [l3] (A) from 33] on mean activity coefficients 7 of NaCl 
as a function of the ionic strength (I) fraction yMgCi 2 of MgCb in NaCl + MgCl 2 mixtures at / = 6 
mol. kg -1 and T = 298.15 K; (B) from ^^] (circles) and 2] (squares) on 7 of NaCl as a function 
of the MgCl 2 molality in NaCl + MgCl 2 mixtures at [NaCl] = 6 mol. kg -1 and T = 298.15 K. 


diminishes exponentially in the hydration shell region Q s h as shown by the profile of 0 PF (r) 
in Fig. 7. The values of a\ , a l 2l and a\ for each activity-concentration curve were obtained 
by first tuning three values of 0(Cf) in Eq. (24) to match three data points In 7 ^) 

with three different concentrations (7 P , j = 1, 2, 3, and then solving the three unknowns 
a\, a 2 , and a\ using three known 0(Cfj) values. For example, for the i = Li + curve in Fig. 
2 A, the selected experimental data points are (y^Ch, In 7 ^-) = (0.315, -0.192), ( 1 , -0.007), 
(1.577, 0.57) and the corresponding tuned 6 *(C P ) are 0.9996, 1.0013, 1.0043. 

The PF model can provide more physical details near the solvated ion (Ca 2+ , for example) 
in a strong electrolyte ([CaCh] = 2 M) such as (1) the dielectric function ?(r) with its varying 
permittivity, (2) variable water density CH 2 o( r )> (3) concentration of counterion C cl -(r), (4) 
electric potential 0 PF (r), and (5) the steric potential 5' trc (r) all shown in Fig. 7. The steric 
potential is small because the configuration of particles (voids between particles) does not 
vary too much from the solvated region to the bulk region. Nevertheless, it has significant 
effect on the variation of mean-held water densities CH 2 o( r ) and hence on the dielectric 
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FIG. 5: Mean activity coefficients of 1:1 electrolyte at various 
results (curves) with experimental data (symbols) compiled in 
coefficients 7 of NaCl in [NaCl] from 0 to 6 mol. kg -1 at T = 
(D) 523.15 (E) 573.15 K. 


temperatures^ Comparison of PF 
{ 13 I from 


291 ] on mean activity 


(A) 298.15 (B) 373.15 (C) 473.15 


function <T(r) in the hydration region. Note that ?(r) is an output, not an input of the 
model. 

The strong electric potential 0 PF (r) in the Born cavity hi (with Rf orn (Cf) = 1.7130 A) 
and the water density CH 2 o( r ) hi the hydration shell Q s h (with ^t 2 + = 5-0769 A) are the 
most important factors allowing the PF results to match the experimental data. The ion and 
shell domains are the crucial region to study ion activities. For example^Fraenkel’s theory is 
entirely based on this region — the so-called smaller-ion shell region [41], The steric energy 
of water molecules modified by the factor Vk+i/vq in Eq. (5) leads to significant changes of 
CH 2 o( r ) and ?(r) profiles in Fig. 7 as compared with those in Fig. 5 in our previous paper 
19]. 
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FIG. 6 : Mean activity coefficients of 2:1 electrolyte at various 
results (curves) with experimental data (symbols) compiled in 
coefficients 7 of MgCl 2 in [MgC^] from 0 to 6 mol. kg ^ 1 at T - 
(D) 523.15 K. 


temperatures. Comparison of PF 


[131 ] from 


3CH32I] on mean activity 


(A) 298.15 (B) 373.15 (C) 423.15 
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Table 3. Values of a[, a 2 , a\ in Eq. (24) 


Fig.# 

i 

a\ a 2 

«3 

Fig.# 

i 

a\ 

a 2 

«3 

2A 

Li+ 

0.99913 0.00069 

0.00009 

3C 

Ca 2+ 

0.99886 

0.00046 

0.00011 

2A 

GI¬ 

0.99893 -0.00008 


3C 

Cl" 

0.99877 

-0.00060 

0.00012 

2B 

LD 

0.99958 -0.00019 

0.00015 

3D 

Ca 2+ 

0.99886 

0.00099 

0.00017 

2B 

Br~ 

0.99822 0.00107 


3D 

Br~ 

0.99920 

-0.00198 

0.00016 

2C 

Na+ 

0.99910 


3E 

Ba 2+ 

0.99844 

0.00011 

0.00010 

2C 

F“ 

0.99933 -0.00029 


3E 

Cl" 

0.99887 

-0.00058 

0.00001 

2D 

Na+ 

0.99927 0.00026 

0.00004 

3F 

Ba 2+ 

0.99851 

0.00054 

0.00008 

2D 

cr 

0.99840 


3F 

Br _ 

0.99926 

-0.00145 

0.00018 

2E 

Na+ 

0.99962 -0.00038 

0.00010 

4A 

Na+ 

1.00581 

-0.00013 


2E 

Br~ 

0.99870 -0.00017 

0.00004 

4B 

Na+ 

1.00527 

0.00042 

0.00019 

2F 

K+ 

0.99934 -0.00120 

0.00007 

5A 

Na+ 

0.9981 


0.0001 

2F 

F“ 

0.99904 0.00013 

0.00004 

5B 

Na+ 

0.9971 

0.0003 

0.0001 

2G 

K+ 

0.99929 -0.00122 

0.00004 

5C 

Na+ 

0.9945 

-0.0007 

0.0001 

2G 

Cl" 

0.99897 -0.00012 

0.00003 

5D 

Na+ 

0.9925 

-0.0028 

0.0001 

2H 

K+ 

0.99931 0.00013 


5E 

Na+ 

0.9870 

-0.0042 

0.0010 

2H 

Br~ 

0.99945 -0.00175 

-0.00006 

6A 

Mg 2+ 

0.9988 

0.0002 

0.0002 

3A 

Mg 2+ 

0.99918 0.00044 

0.00011 

6B 

Mg 2+ 

0.9989 

-0.0004 

0.0003 

3A 

Cl" 

0.99893 -0.00051 

0.00010 

6C 

Mg 2+ 

0.9983 

-0.0014 

0.0005 

3B 

Mg 2+ 

0.99910 0.00063 

0.00015 

6D 

Mg 2+ 

0.9961 

-0.0020 

0.0003 

3B 

Br _ 

0.99888 -0.00065 

0.00018 







Default values: a\ = 1, a 2 = 0, a\ = 0. 


5. CONCLUSION 

A Poisson-Fcrmi model for calculating activity coefficients of aqueous single or mixed 
electrolyte solutions in a large range of concentrations and temperatures has been presented 
and tested by a set of experimental data. The model was shown to well fit experimental 
data with only three adjustable parameters at most for each activity-concentration curve. 
The adjustable parameters correspond to different orders of approximation of the unknown 
Born radius of solvation energy that depends on salt concentrations in a highly complex 



Journal of Chemical Physics 
in the press 


17 


e> 



FIG. 7: Dielectric function ?(r) (denoted by e in the figure), water density CH 2 o( r ) (Ch 2 o)i Cl - 
concentration C ci -(r) ([Cl - ]), electric potential </> PF (r) (0), and steric potential 5 trc (r) (5 trc ) 
profiles near the solvated ion Ca 2+ at [CaCl 2 ] = 2 M, where r is the distance from the center of 
Ca 2+ in Angstrom. 

and nonlinear way. Nevertheless, the values of these parameters have been shown to deviate 
slightly in decimal digits from that of the experimental Born radius in pure water. These 
parameters are physically explained and can be easily verified in future studies for the same 
or different solutions of the present work. The model requires very few parameters because 
it is based on an advanced continuum theory that accounts for steric and correlation effects 
of ions and water with interstitial voids between nonuniform hard spheres. It also deals with 
short- and long-range interactions by partitioning the model domain into the ion, hydration 
shell, and the remaining solvent sub-domains. Numerical methods were also given to show 
how to solve different equations on different sub-domains that describe different physical 
properties of an ion in electrolyte solutions. 
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